################################################################################
#
#       Replication File -- Panel Models
#
#
#       LL, 2.3.2021
#
################################################################################


library(foreign)
library(lfe)
library(lme4)
library(rstanarm)
library(dplyr)
library(arm)
library(tidyr)
library(texreg)
library(brms)


data.border <- read.csv("Data/1_neighbouring_countries.csv")
colnames(data.border)[c(7,30,35,44,51,52)] <- c("Bosnia and Herzegovina", "North Macedonia", "New Zealand", "South Africa", "United Kingdom", "United States")

data1 <- read.csv("Data/2_data_all_merged_std.csv")


# data3 <- data1 %>% 
#   subset(region=="Europe", polity>=6)
data3 <- data1 %>%
  subset(region=="Europe" &  (polity>=6 | is.na(polity)))



standL <- function(x){
  x.out <- (x-mean(x,na.rm=TRUE))/sqrt(var(x,na.rm=TRUE))
  return(x.out)
}


# Standardize some variables
data3$z.restrict_freedom_index <- data3$restrict_freedom_index #DV is not standardized, but I keep new variable name out of simplicity
data3$z.DQ_new <- standL(data3$DQ_new)
data3$z.ConfirmedCases <- standL(data3$ConfirmedCases)
data3$z.ConfirmedCases.perPop <- standL(data3$ConfirmedCases/data3$population)
data3$z.Date <- standL(data3$Date)
data3$z.ConfirmedDeaths <- standL(data3$ConfirmedDeaths)
data3$z.ConfirmedDeaths.perPop <- standL(data3$ConfirmedDeaths/data3$population)
data3$z.hospital_beds <- standL(data3$hospital_beds)
data3$z.gdp_pc <- standL(data3$gdp_pc)
data3$z.RULEOFLAW_new <- standL(data3$RULEOFLAW_new)
data3$z.MUTUCONS_new <- standL(data3$MUTUCONS_new)
data3$z.INDLIB_new <- standL(data3$INDLIB_new)

###########################################################################################
# Generate 10 day intervals

## Taking mean value

data3$month <- substr(data3$Date,5,6)
data3$day <- substr(data3$Date,7,8)
data3 <- data3[-which(is.na(data3$CountryCode)),]

# Creating Time Lag
data3$lag.z.restrict_freedom_index <- NA
for (c in data3$CountryCode){
  pickers <- which(data3$CountryCode==c)
  a <- pickers[7]
  b <- pickers[length(pickers)]
  lag.val <- rep(NA,length(pickers))
  # lag.val[1:10] <- NA
  # for (d in a:b){
  #   within.runner <- which(d==pickers)
  #   lag.val[within.runner] <- mean(data3$z.restrict_freedom_index[pickers][(within.runner-10):(within.runner-1)], na.rm=TRUE)
  # }
  lag.val[1:6] <- NA
  for (d in a:b){
    within.runner <- which(d==pickers)
    lag.val[within.runner] <- data3$z.restrict_freedom_index[pickers][(within.runner-6)]
  }
  data3$lag.z.restrict_freedom_index[pickers] <- lag.val
}

# Creating Time & Spatial Lag
data3$spa.lag.z.restrict_freedom_index <- 0
for (c in data3$CountryCode){
  pickers.for.c <- which(data3$CountryCode==c)
  neighbors.number <- which(data.border[data.border$Neighbouring.countries==data3$CountryName[pickers.for.c[1]],]==1)
  neighbors.names <- colnames(data.border)[neighbors.number]
  pickers.for.neighbors <- which(data3$CountryName %in% neighbors.names)
  for (d in data3$Date[pickers.for.c]){
    day.pick <- which(d==data3$Date)
    aa <- intersect(pickers.for.neighbors,day.pick)
    bb <- intersect(pickers.for.c,day.pick)
    data3$spa.lag.z.restrict_freedom_index[bb] <- mean(data3$lag.z.restrict_freedom_index[aa], na.rm = TRUE)
  }
}






data4 <- data3 %>%  
  mutate(m.period =  case_when(
    day > 0 & day<11 ~ 1,
    day > 10 & day<21 ~ 2,
    day > 20  ~ 3))

data4 <- data4 %>%
  mutate(date.new = as.numeric(paste0(as.character(month),as.character(m.period))))


# data4.ss <- data4[,c(1,2,154:162)] 
# 
# 
# data5 <- data4.ss %>%
#   group_by(m.period,CountryName,month) %>%
#   summarise(z.restrict_freedom_index = mean(z.restrict_freedom_index),
#                z.DQ_new = mean(z.DQ_new, na.rm = TRUE),
#                z.ConfirmedCases = mean(z.ConfirmedCases, na.rm = TRUE),
#             z.ConfirmedDeaths = mean(z.ConfirmedDeaths, na.rm = TRUE),
#             date.new = mean(date.new, na.rm = TRUE))
#                
# data5$z.date.new <- standL(data5$date.new)



##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##
### Daily Models


### Daily Models
mod9 <- lmer(z.restrict_freedom_index ~ z.INDLIB_new +  z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod9)


mod6 <- lmer(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+ 
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod6)



mod7 <- lmer(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod7)


mod8 <- lmer(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod8)



screenreg(list(mod7,mod9,mod8,mod6), 
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
          reorder.coef  = c(2,11,12,13,3,4,8,5,6,7,9,10,1),
          custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"))

htmlreg(list(mod7,mod9,mod8,mod6), 
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
        reorder.coef  = c(2,11,12,13,3,4,8,5,6,7,9,10,1),
        custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/Table4_main.html")



### Country-wise drop 

countries <- unique(attributes(mod9)$frame$CountryName)

coef.catcher <- matrix(NA,length(countries),4)
var.catcher <- matrix(NA,length(countries),4)
for (i in countries){
  ### Daily Models
  mod9 <- lmer(z.restrict_freedom_index ~ z.INDLIB_new +  z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
                 lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
                 z.Date + I(z.Date*z.Date)+ 
                 I(z.Date*z.Date*z.Date) + 
                 (1|CountryName),
               data = data3[data3$CountryName!=i,])
  
  
  
  mod6 <- lmer(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+ 
                 lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
                 z.Date + I(z.Date*z.Date)+ 
                 I(z.Date*z.Date*z.Date) + 
                 (1|CountryName),
               data = data3[data3$CountryName!=i,])
  
  
  
  
  mod7 <- lmer(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
                 lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
                 z.Date + I(z.Date*z.Date)+ 
                 I(z.Date*z.Date*z.Date) + 
                 (1|CountryName),
               data = data3[data3$CountryName!=i,])
  
  
  
  mod8 <- lmer(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
                 z.Date + I(z.Date*z.Date)+ 
                 I(z.Date*z.Date*z.Date) + 
                 (1|CountryName),
               data = data3[data3$CountryName!=i,])
  
  coef.catcher[which(i==countries),1] <-  fixef(mod9)[2]
  coef.catcher[which(i==countries),2] <-  fixef(mod6)[2]
  coef.catcher[which(i==countries),3] <-  fixef(mod7)[2]
  coef.catcher[which(i==countries),4] <-  fixef(mod8)[2]
  
  var.catcher[which(i==countries),1] <-  vcov(mod9)[2,2]
  var.catcher[which(i==countries),2] <-  vcov(mod6)[2,2]
  var.catcher[which(i==countries),3] <-  vcov(mod7)[2,2]
  var.catcher[which(i==countries),4] <-  vcov(mod8)[2,2]
}


z.catcher <- coef.catcher/sqrt(var.catcher)

pdf(file = "Out/FigureA4_app.pdf", 
    width = 9, 
    height = 6)
par(oma=c(6,1,0,0))
plot(0, xaxt = 'n',  bty = 'n', pch = '', ylab = '', xlab = '', 
     xlim=c(0,33), ylim=c(-3,.5))
for (i in 1:length(countries)){
  points(i,coef.catcher[i,1], pch=19)
  segments(i,coef.catcher[i,1]-1.96*sqrt(var.catcher)[i,1],i,coef.catcher[i,1]+1.96*sqrt(var.catcher)[i,1], lwd=0.5)
  segments(i,coef.catcher[i,1]-1.64*sqrt(var.catcher)[i,1],i,coef.catcher[i,1]+1.64*sqrt(var.catcher)[i,1], lwd=1.5)
}
abline(h=0, lwd=2)
axis(side =1, at = c(1:32),las = 3,labels = countries)
dev.off()




######## Without GDP


### Daily Models
mod9 <- lmer(z.restrict_freedom_index ~ z.INDLIB_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + 
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod9)


mod6 <- lmer(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+ 
               lag.z.restrict_freedom_index+ z.Date + 
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod6)



mod7 <- lmer(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + 
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod7)


mod8 <- lmer(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + 
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod8)



screenreg(list(mod7,mod9,mod8,mod6), 
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
          reorder.coef  = c(2,10,11,12,3,4,5,6,7,8,9,1),
          custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time",  "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"))

htmlreg(list(mod7,mod9,mod8,mod6), 
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
        reorder.coef  = c(2,10,11,12,3,4,5,6,7,8,9,1),
        custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time",  "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA6_app.html")


# Sub Sample

data_vdem <- data3 %>% 
  subset(CountryName!="Serbia" & (CountryName!="Hungary")& (CountryName!="Albania"))


### Daily Models
mod9 <- lmer(z.restrict_freedom_index ~ z.INDLIB_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data_vdem)
summary(mod9)


mod6 <- lmer(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+ 
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data_vdem)
summary(mod6)



mod7 <- lmer(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data_vdem)
summary(mod7)


mod8 <- lmer(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data_vdem)
summary(mod8)



screenreg(list(mod7,mod9,mod8,mod6), 
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
          reorder.coef  = c(2,11,12,13,3,4,8,5,6,7,9,10,1),
          custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"))

htmlreg(list(mod7,mod9,mod8,mod6), 
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
        reorder.coef  = c(2,11,12,13,3,4,8,5,6,7,9,10,1),
        custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA9_app.html")


### with backsliding

data3$libdem10 <- -data3$libdem10

mod9 <- lmer(z.restrict_freedom_index ~ z.INDLIB_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+ 
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + libdem10 + 
               (1|CountryName),
             data = data3)
summary(mod9)


mod6 <- lmer(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+ 
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + libdem10 + 
               (1|CountryName),
             data = data3)
summary(mod6)



mod7 <- lmer(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + libdem10 + 
               (1|CountryName),
             data = data3)
summary(mod7)


mod8 <- lmer(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + libdem10 + 
               (1|CountryName),
             data = data3)
summary(mod8)



screenreg(list(mod7,mod9,mod8,mod6), 
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
          reorder.coef  = c(2,12,13,14,3,4,8,5,6,7,9,10,11,1),
          custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", 
                                "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Backsliding", 
                                "Individual Liberties", "Rule of Law", "Mutual Constraints"))

htmlreg(list(mod7,mod9,mod8,mod6), 
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
        reorder.coef  = c(2,12,13,14,3,4,8,5,6,7,9,10,11,1),
        custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", 
                              "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Backsliding", 
                              "Individual Liberties", "Rule of Law", "Mutual Constraints"),   
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA12_app.html")





### w/o democracy variables


mod8 <- lmer(z.restrict_freedom_index ~  z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod8)



screenreg(list(mod8), 
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Baseline Model"),#, "Model 3.2", "Model 3.3", "Model 3.4"),
          reorder.coef  = c(2,3,7,4,5,6,8,9,1),
          custom.coef.names = c("Intercept","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3"))

htmlreg(list(mod8), 
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Baseline Model"),
        reorder.coef  = c(2,3,7,4,5,6,8,9,1),
        custom.coef.names = c("Intercept","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA14_app.html")








### Interaction models for diffusion
# Note: The lme4 functions lead sometimes to convergence warnings. To rule out problem we re-estimate the models
# in the Bayesian framework and ensure that point estimates are similar.
moda9 <- lmer(z.restrict_freedom_index ~ z.INDLIB_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+  
                spa.lag.z.restrict_freedom_index*z.INDLIB_new+ z.Date + z.gdp_pc+
                z.Date + I(z.Date*z.Date)+ 
                I(z.Date*z.Date*z.Date) + 
                (1+spa.lag.z.restrict_freedom_index|CountryName),
              data = data3)

# to check since convergence warnings
# modb9 <- brm(z.restrict_freedom_index ~ z.INDLIB_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+  
#                spa.lag.z.restrict_freedom_index*z.INDLIB_new+ z.Date + z.gdp_pc+
#                z.Date + I(z.Date*z.Date)+ 
#                I(z.Date*z.Date*z.Date) + 
#                (1+spa.lag.z.restrict_freedom_index|CountryName),
#              data = data3)
# 
# summary(moda9)
# summary(modb9)



moda6 <- lmer(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+ 
                spa.lag.z.restrict_freedom_index*z.MUTUCONS_new+ z.Date + z.gdp_pc+
                z.Date + I(z.Date*z.Date)+ 
                I(z.Date*z.Date*z.Date) + 
                (1+spa.lag.z.restrict_freedom_index|CountryName),
              data = data3)

# modb6 <- brm(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+ 
#                spa.lag.z.restrict_freedom_index*z.MUTUCONS_new+ z.Date + z.gdp_pc+
#                z.Date + I(z.Date*z.Date)+ 
#                I(z.Date*z.Date*z.Date) + 
#                (1+spa.lag.z.restrict_freedom_index|CountryName),
#              data = data3)
# 
# summary(moda6)
# summary(modb6)


moda7 <- lmer(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+  
                spa.lag.z.restrict_freedom_index*z.DQ_new+ z.Date + z.gdp_pc+
                z.Date + I(z.Date*z.Date)+ 
                I(z.Date*z.Date*z.Date) + 
                (1+spa.lag.z.restrict_freedom_index|CountryName),
              data = data3)

# modb7 <- brm(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+  
#                spa.lag.z.restrict_freedom_index*z.DQ_new+ z.Date + z.gdp_pc+
#                z.Date + I(z.Date*z.Date)+ 
#                I(z.Date*z.Date*z.Date) + 
#                (1+spa.lag.z.restrict_freedom_index|CountryName),
#              data = data3)
# summary(moda7)
# summary(modb7)




moda8 <- lmer(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+  
                spa.lag.z.restrict_freedom_index*z.RULEOFLAW_new+ z.Date + z.gdp_pc+
                z.Date + I(z.Date*z.Date)+ 
                I(z.Date*z.Date*z.Date) + 
                (1+spa.lag.z.restrict_freedom_index|CountryName),
              data = data3)

# modb8 <- brm(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + lag.z.restrict_freedom_index+  
#                spa.lag.z.restrict_freedom_index*z.RULEOFLAW_new+ z.Date + z.gdp_pc+
#                z.Date + I(z.Date*z.Date)+ 
#                I(z.Date*z.Date*z.Date) + 
#                (1+spa.lag.z.restrict_freedom_index|CountryName),
#              data = data3)
# summary(moda8)
# summary(modb8)



screenreg(list(moda7,moda9,moda8,moda6), 
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
          reorder.coef  = c(2,11,12,13,14,15,16,17,6,3,4,8,5,7,9,10,1),
          custom.coef.names = c("Intercept", "Quality of Democracy", "Confirmed Deaths (rel.)", "# of Hospital Beds", "Lagged Outcome (Y_[t-7])",
                                "Spatial Lag (Neighbors)", "Time", "GDP pc", "Time^2", "Time^3", "Quality of Democracy X Sp Lag",
                                "Individual Liberties", "Individual Liberties X Sp Lag", "Rule of Law ","Rule of Law X Sp Lag",
                                "Mutual Constraints", "Mutual Constraints X Sp Lag"))


htmlreg(list(moda7,moda9,moda8,moda6), 
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
        reorder.coef  = c(2,11,12,13,14,15,16,17,6,3,4,8,5,7,9,10,1),
        custom.coef.names = c("Intercept", "Quality of Democracy", "Confirmed Deaths (rel.)", "# of Hospital Beds", "Lagged Outcome (Y_[t-7])",
                              "Spatial Lag (Neighbors)", "Time", "GDP pc", "Time^2", "Time^3", "Quality of Democracy X Sp Lag",
                              "Individual Liberties", "Individual Liberties X Sp Lag", "Rule of Law ","Rule of Law X Sp Lag",
                              "Mutual Constraints", "Mutual Constraints X Sp Lag"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA15_app.html")











### Alternative data 2017 DB data


data1 <- read.csv("Data/2_data_all_merged_robust_std.csv")


# data3 <- data1 %>% 
#   subset(region=="Europe", polity>=6)
data3 <- data1 %>%
  subset(region=="Europe" &  (polity>=6 | is.na(polity)))



standL <- function(x){
  x.out <- (x-mean(x,na.rm=TRUE))/sqrt(var(x,na.rm=TRUE))
  return(x.out)
}


# Standardize some variables
data3$z.restrict_freedom_index <- data3$restrict_freedom_index #DV is not standardized, but I keep new variable name out of simplicity
data3$z.DQ_new <- standL(data3$DQ_new)
data3$z.ConfirmedCases <- standL(data3$ConfirmedCases)
data3$z.ConfirmedCases.perPop <- standL(data3$ConfirmedCases/data3$population)
data3$z.Date <- standL(data3$Date)
data3$z.ConfirmedDeaths <- standL(data3$ConfirmedDeaths)
data3$z.ConfirmedDeaths.perPop <- standL(data3$ConfirmedDeaths/data3$population)
data3$z.hospital_beds <- standL(data3$hospital_beds)
data3$z.gdp_pc <- standL(data3$gdp_pc)
data3$z.RULEOFLAW_new <- standL(data3$RULEOFLAW_new)
data3$z.MUTUCONS_new <- standL(data3$MUTUCONS_new)
data3$z.INDLIB_new <- standL(data3$INDLIB_new)

###########################################################################################
# Generate 10 day intervals

## Taking mean value

data3$month <- substr(data3$Date,5,6)
data3$day <- substr(data3$Date,7,8)
data3 <- data3[-which(is.na(data3$CountryCode)),]

# Creating Time Lag
data3$lag.z.restrict_freedom_index <- NA
for (c in data3$CountryCode){
  pickers <- which(data3$CountryCode==c)
  a <- pickers[7]
  b <- pickers[length(pickers)]
  lag.val <- rep(NA,length(pickers))
  lag.val[1:6] <- NA
  for (d in a:b){
    within.runner <- which(d==pickers)
    lag.val[within.runner] <- data3$z.restrict_freedom_index[pickers][(within.runner-6)]
  }
  data3$lag.z.restrict_freedom_index[pickers] <- lag.val
}

# Creating Time & Spatial Lag
data3$spa.lag.z.restrict_freedom_index <- 0
for (c in data3$CountryCode){
  pickers.for.c <- which(data3$CountryCode==c)
  neighbors.number <- which(data.border[data.border$Neighbouring.countries==data3$CountryName[pickers.for.c[1]],]==1)
  neighbors.names <- colnames(data.border)[neighbors.number]
  pickers.for.neighbors <- which(data3$CountryName %in% neighbors.names)
  for (d in data3$Date[pickers.for.c]){
    day.pick <- which(d==data3$Date)
    aa <- intersect(pickers.for.neighbors,day.pick)
    bb <- intersect(pickers.for.c,day.pick)
    data3$spa.lag.z.restrict_freedom_index[bb] <- mean(data3$lag.z.restrict_freedom_index[aa], na.rm = TRUE)
  }
}






data4 <- data3 %>%  
  mutate(m.period =  case_when(
    day > 0 & day<11 ~ 1,
    day > 10 & day<21 ~ 2,
    day > 20  ~ 3))

data4 <- data4 %>%
  mutate(date.new = as.numeric(paste0(as.character(month),as.character(m.period))))



##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##    ##
### Daily Models
mod9 <- lmer(z.restrict_freedom_index ~ z.INDLIB_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod9)


mod6 <- lmer(z.restrict_freedom_index ~ z.MUTUCONS_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+ 
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod6)



mod7 <- lmer(z.restrict_freedom_index ~ z.DQ_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  
               lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod7)


mod8 <- lmer(z.restrict_freedom_index ~ z.RULEOFLAW_new + z.ConfirmedDeaths.perPop+ z.hospital_beds + spa.lag.z.restrict_freedom_index+  lag.z.restrict_freedom_index+ z.Date + z.gdp_pc+
               z.Date + I(z.Date*z.Date)+ 
               I(z.Date*z.Date*z.Date) + 
               (1|CountryName),
             data = data3)
summary(mod8)



screenreg(list(mod7,mod9,mod8,mod6), 
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
          reorder.coef  = c(2,11,12,13,3,4,8,5,6,7,9,10,1),
          custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"))

htmlreg(list(mod7,mod9,mod8,mod6), 
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 3.1", "Model 3.2", "Model 3.3", "Model 3.4"),
        reorder.coef  = c(2,11,12,13,3,4,8,5,6,7,9,10,1),
        custom.coef.names = c("Intercept", "Quality of Democracy","Confirmed Deaths (rel.)", "# of Hospital Beds", "Spatial Lag (Neighbors)", "Lagged Outcome (Y_[t-7])","Time", "GDP pc", "Time^2", "Time^3", "Individual Liberties", "Rule of Law", "Mutual Constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA3_app.html")



